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Abstract. In this paper we construct Discontinuous Galerkin approximations of the Stokes 
problem where the velocity field is H (div, f2)-conforming. This implies that the velocity 
solution is divergence-free in the whole domain. This property can be exploited to design a 
simple and effective preconditioner for the final linear system. 



1. Introduction 

In this paper we present a preconditioning strategy for a family of discontinuous Galerkin 
H (div; f2)-conforming discretizations of the Stokes problem. The methods considered here 
were introduced in [TU] where the authors showed that the approximate velocity field is truly 
divergence-free. These same methods were also used in [H]. 

In general, the numerical discretization of the Stokes problem produces algebraic linear 
systems of equations of the saddle-point type. Solving such algebraic linear systems has been 
the subject of considerable attention from various communities and nowadays at present 
several different approaches can be used to solve them efficiently (see [UJ and references 
cited therein). One widely used approach is to use a block diagonal preconditioner with 
two blocks: one containing the inverse or a preconditioner of the stiffness matrix of a vector 
Poisson discretization, and one containing the inverse of a lumped mass matrix for the 
pressure. The resulting preconditioned system is then solved by means of the MINRES 
(Minimal Residual). For this type of preconditioner, a sound mathematical theory can be 
developed that guarantees uniform convergence (in the asymptotic limit) with respect to 
the mesh size of the discretization. However, it is also generally observed that with such a 
preconditioner the convergence of the MINRES can be quite slow. 

One possible reason for the slow convergence is related to the treatment of the divergence- 
free condition on the velocity field. As is well known, the question of how to satisfy such a 
condition at the discrete level is one of the most intricate in Finite Elements. 

Herein, we consider the family of discontinuous Galerkin H (div; f2)-conforming methods 

introduced in [TO]. When appropriate finite element spaces are selected, the resulting ap- 
proximation to the velocity field turns out to be exactly divergence-free. This property can 
be fully exploited in designing and constructing efficient preconditioners. Indeed, by using it, 
we can reduce the solution of the Stokes problem to the solution of a "second-order" problem 
in the space cuyXHq (fl). In this way, the incompressibility constraint is incorporated into 
the linear system to be solved, such that the linear system is no longer in saddle-point form. 
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We propose then a preconditioner for the solution of the corresponding problem in curl Hq(Q). 
This is done by means of the fictitious space [HU [TT] (or auxiliary space (2TJJ, [18]) frame- 
work, and the proposed preconditioner is the solution of several elliptic problems: a vector 
Laplacian discretized with (div; fi)-conforming methods, and another elliptic problem 

discretized with an if ^conforming method. The solution of such systems can then be effi- 
ciently computed with classical approaches, for instance the Geometric Multigrid (GMG) or 
Algebraic Multigrid (AMG) methods. 

Throughout the paper, we use the standard notation for Sobolev spaces [1] . For a bounded 
domain D C M 2 , we denote by H m (D) the L 2 -Sobolev space of order m > and by || • || mi £> 
and | ■ \ m ,D the usual Sobolev norm and seminorm, respectively. For m = 0, we write L 2 (D) 
instead of H°(D). For a general summability index p, we also denote by W m ' p (D) the usual 
L p -Sobolev spaces of order m > with norm || ■ || mjP) D and seminorm | ■ | m)Pl i> By conven- 
tion, we use boldface type for the vector-valued analogues: H m (D) = [H m (D)} 2 , likewise, 
we use boldface italics for the symmetric-tensor-valued analogues: g H m (D) := [H m {D)] 2 ^. 
H m (D)/M. denotes the quotient space consisting of equivalence classes of elements of H m (D) 
that differ by a constant; for m = the quotient space is denoted by L 2 (D)fR. We indicate 
by Ll(D) the space of the L 2 (D) functions with zero average over D (which is obviously iso- 
morphic to L 2 (D)/M.). We use (■ ,-)d to denote the inner product in the spaces L 2 (D), L 2 (D), 
and C 2 (D). 

Let Q C M 2 be a polygonal domain with boundary T = dQ, and let / G [L 2 (Q)] 2 . The 
Stokes equations for a viscous incompressible fluid are 



(1.1) 



-dW(2ue(u)) + Vp = f in Q 
div u = in Q 



where, with the usual notation, u is the velocity field, p the pressure, v the viscosity of 
the fluid, and e{u) G [£ 2 (^)]sym the symmetric (linearized) strain rate tensor defined by 
s(u) = !(Vw + (Vtt) T ). Moreover, we assume no-slip boundary conditions 

(1.2) un = onT, 



together with the natural condition associated with (1.1) on the tangential component of 
the normal stresses 

(1.3) ((2i/e(«)-pl)-n)-t = on T, 



where I is the identity tensor. Note that as (I • n) • t = then (1.3) is reduced to 

(1.4) (e(u) • n) • t = on T. 

When the space H Qn (Q) = {v G H (Q) : v ■ n = on T } is introduced, the variational 
formulation of the Stokes problem reads: Find (u,p) G H Qn (Q,) x L 2 (fi)/IR as the solution 



(1.5) 
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of: 

;(u, v) + b(v,p) = / f-vdx Vw G ifj )ri (ft) 
Jn 

b(u,q) = Vg 6 L 2 (n)/1 

where the bilinear forms are defined by 

a(u,v) := 2u / e(u) : £(v) dx, Vu, v G .tfjjj >n (ft), 

6( v>9 ) : = - ! gdivvdx Vg G L 2 (ft)/M ,Vu G ffj in (fi). 



The existence and uniqueness of the solution (u,p) of (1.5) are very well known (see for 
instance [EH H21 D3])- 

2. Abstract setting and basic notations 

Let Th be a shape-regular family of partitions of ft into triangles T. We denote by the 
diameter of T, and we set h = maxyg^ We also assume that the decomposition Th is 
conforming in the sense that it does not contain any hanging nodes. 

We denote by Eh the set of all edges and by E% and E% the collection of all interior and 
boundary edges, respectively. 

For s > 1, we define 

H s (T h ) = {4>e L 2 (ft) , such that 4>\ T G H S (T), VT G T h } , 

and their analogues H s (Th) and 'H s (Th), respectively. For scalar, vector- valued, and tensor 
functions, we use (• , -)<j- h to denote the L 2 (7h)-inner product and (■ , -)g h to denote the L 2 (£ h )- 
inner product elementwise. 

The vector functions are represented column-wise. We recall the definitions of the following 
operators acting on vectors v G if 1 (ft) and on scalar functions G H 1 ^) as 

dv 1 , dv 2 dv 1 r 1 > i < 



(2.1) divv = > - — , curlu = — — , cur\(f) = V ± 4 

dxi ax\ 0x2 

i=i 

And, we recall the definitions of the spaces to be used herein: 

H(div; ft) := {v G X 2 (ft) : divv G L 2 (ft) }, 
if (curl; ft) := {v G £ 2 (ft) : curl v G L 2 (ft) }. 



90 dcj) 
8x2 ' <9xi 



i?o,n(div; ft) := G ff(div; ft) : v ■ n = on T }, 
ff ,n(div°;ft) := G if ,n(div;ft) : divv = in ft}. 
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The above spaces are Hilbert spaces with the norms 

IMI;j(div,n) : = \\ v \\o,n + l|div«||o, n Vu G ff(div; Q) , 
IMIff(curi,n) : = Hlo,n + Hcurlwlljn G if (curl; fi). 

Remark 2.1. zs worth noting that if we restrict our analysis to vectors u and v in 
i? 1 (fi) PI Ho !n (div ; f2) i/ien problem (1.5) becomes: Find u G i? 1 (fi) PI Ho >n (div°; ft) as the 



solution of: 

(2.2) a{u,v) = [ f-vdx Vv G if^Q) n H 0>n (div°; Q). 

Jn 

As is usual in the DG approach, we now define some trace operators. Let e G £% be 
an internal edge of Tu shared by two elements T 1 and T 2 , and let n 1 (n 2 ) denote the unit 
normal on e pointing outwards from T 1 (T 2 ). For a scalar function ip G -?/ 1 (7/ l ), a vector 
field r G H^Th), 

or a tensor field r G / H 1 (Th) we define the average operator in the usual 
way (see for instance |3]), that is, on internal edges 

M = + V% {v} = \{v X + v% {t} = ^(t 1 + t 2 ). 

However, on a boundary edge, we take {ip}, {v}, and {r} as the trace of ip, v, and r 
Respectively, on that edge. 

For a scalar function ip G if 1 (7/ l ), the jump operator is defined as 

[ ip ] = y^n 1 + (/9 2 n 2 on e G £%, and [ v 3 ] = V 911 on e G £f 

(where obviously n is the outward unit normal), so that the jump of a scalar function is a 
vector in the normal direction. 

For a vector field v G if 1 (7^), following, for example, jl], the jump is the symmetric 
matrix-valued function given on e by 

{vj = v 1 n 1 + v 2 n 2 on e G ££, and [«] = v n on e G 

where vQn = (vn T -\-nv T )/2 is the symmetric part of the tensor product of v and n. Hence, 
the jump of a vector-valued function is a symmetric tensor. 

If we denote by n-r the outward unit normal to dT, it is easy to check that 

(2.3) [ v-n T qds=J2[{ v }-l ( l} ds V v e H(div; Q) , Vg G H\%). 

TeT h JdT ee£ h Je 

Also for r G and for all t> G if 1 (7^), we have 

(2.4) £ / (T.n r ).«d a =X;/{r}:Md a . 
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2.1. Discrete Spaces: General framework. We present three choices for each of the 
finite element spaces Vh and Qh to approximate velocity and pressure, respectively. For 
each choice, we also need an additional space Mh made of piecewise polynomial scalars to be 
used as a sort of potentials. We will explain the reason for doing this and the way in which 
to do this later on. Note, too, that we will use this space more heavily in the construction 
of our preconditioner. The different choices for the spaces Vh, Qh, and Mh rely on different 
choices of the local polynomial spaces 1Z{T), S(T), and A4(T), respectively, made for each 
element T. Specifically, we have 

(2.5) V h := {v E H{div; Q) : v\ T E ll(T) V T E T h , v ■ n = on T} , 

(2.6) Q h := {q E L 2 (n)/R : q\ T E S(T) VT E %} , 
and 

(2.7) M h := {<p E H^O) : <p\ T E M(T) V T E %} . 

The three spaces Vh, Qh, and Mh will always be related by this exact sequences: 

(2.8) 0^M h ^n^Q^O. 



It is also necessary for each operator in (2.8) to have a continuous right inverse whose norm 



is uniformly bounded in h. For instance, it is necessary that 

(2.9) 3/3 > s.t. \/h,\/q E Q h Bv E V h with: divv = q and ||«||o,n < dklkfi- 

Remark 2.2. In all our examples, the pair (Vh, Qh) is among the classical (and very old) 
finite element spaces specially tailored for the approximation of the Poisson equation in mixed 



form. In particular, properties (2.8) and (2.9) always hold. 



2.2. Examples. We now present three examples of finite element spaces that can be used 
in the above framework. For each example, we specify the corresponding polynomial spaces 
used on each element and describe the corresponding sets of degrees of freedom. We restrict 
our analysis to the case of triangles; the case of parallelograms can also be considered when 
corresponding changes are made(see [H]). 

Let us first fix the notation concerning the spaces of polynomials. For m > 0, we denote 
by P m (T) the space of polynomials defined on T of degree of at most m; the corresponding 
vector space is denoted by P m (T) = (P m (T)) 2 . A polynomial of degree m > 3 that vanishes 
throughout dT (hence it belongs to Hq(T)) is called a bubble (or an H-bubble) of degree m 
over T . The space of bubbles of degree m over T is denoted by HB m (T). We denote by 
P™ om (r) the space of homogeneous polynomials of degree m, and we denote by x- 1 the vector 
(-x 2 ,x 1 ). 

For m > 2, 

(2.10) P+(T) := P m (T) + HB m+l {T) 
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And, for m > 1, we set 

(2.11) BDM m (T) := P m (T), RT m (T) := P m (T) © xP™ m (T) 
and 

(2.12) TR m (T) := P m (T) © x^P&JT). 

We also consider some generalized bubbles: a vector-valued polynomial of degree m > 2 that 
belongs to i?o,n(div, T) (hence whose normal component vanishes throughout dT) is called 
a D-bubble of degree m over T. The space of D-bubbles of degree m over T is denoted by 
DB m {T). 

All the spaces used herein are well known and widely used. They are usually referred to as 
Brezzi-Douglas-Marini, Raviart- Thomas, and Rotated Raviart- Thomas spaces, respectively. 
The first example follows. 

1. Raviart- Thomas For k > 1, we take in each T, S(T) = P fe (T), and 1Z(T) := RT fe (T). 
The degrees of freedom in RTfc(T) are 



(2.13) 



J^u-n e qds Ve G dT, Vg G P fe (e), 
/ u-pdx VpGP^(T). 



As Qh is made of discontinuous piecewise polynomials, here and in the following examples 
the degrees of freedom in S (T) can be taken in an almost arbitrary way. The corresponding 
pair of spaces (Vh, Qh) gives the classical Raviart-Thomas finite element approximation for 
second-order elliptic equations in mixed form, as introduced in [19|. It is well known and 
easy to check that the pair (Vh, Qh) satisfies 

(2.14) div(V h ) = Q h 



and that the property (2.9) is verified. We then take M.(T) := P fc+1 (T), and note that 

(2.15) curl(A4) C V h - 

We also note that the operator curl has a continuous right inverse uniformly bounded from 
V h H i3" ,n(div°, CI) to A4; that is, 

3 C > such that V7i, Vv h G Vh H H 0)n (div°, Cl) 3^ G M h , such that 

(2.16) cur\(p = v h and |M|i,n < C \\v h \\o,u- 

2. Brezzi-Douglas-Marini: For k > 1, we take «S(T) = P fc_1 (T), and 1Z(T) = BDM fc (T). 
The degrees of freedom for BDMfc(T) are (see [5]): 



(2.17) 



^u-n e qds Ve G dT, M q G P fc (e); 

J u-vdx Vw G TR fc _ 2 (T) fe>2. 
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The resulting finite element pair (Vh, Qh) is also commonly used for the approximation of 
second-order elliptic equations in mixed form introduced in [8]. Also in this case it has been 
established that the pair (Vh, Qh) verifies the properties of (2.14) and (2.9). We then take 
M{T):= F k+1 (T), and note that dzTHb and flOjb are also satisfied. 



3. Brezzi-Douglas-Fortin-Marini: For A; > 1, we take <S(T) = F k (T) and H(T) = BDFM fe (T), 
which can be written as BDFMfc = BDMj,(T) + DBk + i(T). The degrees of freedom for 
BDFMfc(T), though similar to the previous ones, are given here: 



J^u-n e qds We G dT, Wq G P fc (e); 
u-vdx VoGTRn(T). 



(2.18) 



The resulting finite element pair (V^, Qh) gives the triangular analogue of the element 
BDFMfc introduced in [7] for the approximation of second-order elliptic equations in mixed 



form. It is easy to check that the pair (Vh, Qh) verifies (2.14) and (2.9). We then take 



M(T) : = P£ +1 (T), and note that ( |2T5| ) and ( gig ) hold. 



The three choices above are quite similar to each other, and the best choice among them 
generally depends on the problem and the way in which the discrete solution is to be used. 
We also use basic approximation properties: for instance, we recall that a constant C exists 
such that for all T G Th and for all v, e.s. in H S (T), an interpolant v 1 G 7t(T) exists such 
that 



(2.19) 



11^ — v \\o,t + ^t|^ 7 |i,t < Ch s T \v\ Sj Ti s < k + 1. 



3. The discontinuous Galerkin H(dtv; ^)-conforming method 

To introduce our DG-approximation, we start by defining, for any u, v G H 2 {T~h) and any 
p, q G L 2 (Q)fR, the bilinear forms 



A h (u,v) = 2u [(e(u) : e(v)) Th - <{*(«)} : [v]) E o - ([«] : {e(v)}) S o] 



(3.1) 



2u 



+ 2v 



(e(u)n, (v ■ n)n) g 8 + ({u ■ n)n, e(v)n) e 



J2 ah e' /[«] = [^Ids+^a/i, 1 [{u-n){vn)ds 



5 h («,g) = -(g,divi;) rh V« G H 2 {T h )^q G L 2 (i2)/I 



where as usual a is the penalty parameter that we assume to be positive and large enough. 
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It is easy to check that the solution (u,p) of ( |1.5 ) verifies: 



(3.2) 



A h (u, v) + B h (v, p) = (f, v) 
B h {u,q) = 



Wv e H 2 {%) 

Vg G L 2 (n)/R. 



For a general DG approximation, we now replace the spaces H 2 {T~h) and L 2 (f2)/IR with the 
discrete ones Xh and Qh, respectively. Following [TD], we choose for (X^, Qh) one of the pairs 
(Vh, Qh) of the previous examples in order to get a global divergence-free approximation. 

More generally, we can choose a pair (V^, Qh) in order to find a third space Mh in such 
a way that (2.8), (2.14), (2.9), (2.15), and (2.16) are satisfied. This set of assumptions will 
come out several times in the sequel and, therefore, it is helpful to give it a special name. 

Definition 3.1. In the above setting, we say that the three spaces (Vh, Qh,Nh) satisfy As- 
sumption HO if ( pT8j ), ( |2l4| ), fl2~9j ), ( [215] and fl2~l6j ) are satisfied. 

We note that, according to the definition of V^, the normal component of any v G is 
continuous on the internal edges and vanishes on the boundary edges. Therefore, by splitting 
a vector v G Vh into its tangential and normal components v n and v t 



(3.3) 
we have 

(3.4) 

implying that 
(3.5) 



v, 



Ve G £ h 



[v ■ n n, 



(v ■ t)t = v — v r 



rds = Vr G ^H l (Th), 



Vee£ h Jlv]:rds = J[v t ]:rds VreU x {T h ). 



The resulting approximation to (1.5), therefore, becomes: Find (u h ,Ph) in Vh x Qh such 
that 



(3.6) 
where 

(3.7) 



a h (u h , v) + b(v, ph) = (/, v) 
b(u h , q) = 

:= 2i/ [(e(u) : e(v)) Th - ({e(u)} 



Vv G V h 

VgG Q fc , 



:«*]:{e(«)}>fi-] 



b(v, q) := — (q, divv) 



ut]:[v t ]ds \/u,veV h , 
VveV h , VqeQh. 



Consistency The consistency of the formulation (3.6 ) can be checked by means of the usual 
DG-machinery. In this case, it is sufficient to compare (3.1) and (3.7) and to observe that if 
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(u,p) is the solution of (1.5), then 

A h (u,v h ) = a h (u,v h ), B h (v h ,p) = b(v h ,p), Vv h eV h C H 0jn (dbr,Q), 

Further, it is evident that, Bh(u,qh) = b(u,qh) for all qh G Qh- Hence, as (u,p) verifies 
(3.2), it also verifies (3.6); that is, 

a h (u,v) + b(v,p) = (/,«) \Jv G V h 

b(u,q) = VqeQh- 



(3.8) 



Thus, consistency is proved. 

To prove the existence and uniqueness of the solution of (3.6) and to obtain the optimal 
error bounds, we need to define suitable norms. We define the following semi-norms 



0,T5 



IK 



I>.-i 



v 



I0,e> 



TeT h 



and norms 



(3.9) 



v 



DG 



\\v\\ 



v G H\T h \ 
v\\dg+ J2 2ph T\ £ ( v ">\lT veH 2 {T h ). 



TeT h 



And, we remark that the seminorms defined in (3.9) are actually norms with the addi- 
tional requirement that v G -Ho,n (div; fl). We also observe that when restricted to discrete 
functions v G Vh, the || • H^G-norm and the ||| • ||| are equivalent (using inverse inequality). 
Continuity can easily be shown for both bilinear forms: 

\a h (u,v)\ < \\u\\ 1 v 1 Vv, v G H 2 (T h ), 

\b(v,q)\ < ||«||i jft ||g|| 0> n V« G H\%), q G L 2 (n)/R . 

Following [9], the existence and uniqueness of the approximate solution and optimal error 
bounds are guaranteed if the following two conditions are satisfied: 

(HI): coercivity: 3 7 > independent of the mesh size h such that 
(3.10) a h (v,v) > l\\v\\ 2 DG VveV h . 

(H2): inf-sup condition: 3 > independent of the mesh size h such that 

(divv,q h ) n 



(3.11) 



sup 
vev h 



\ v \\dg 



> P\\qh\\o,n Vqh G Q h . 



Condition (H2) is a consequence of the inf-sup condition that holds for the continuous 
problem ( |1.5 ): 

3/3 > s.t. V/i, \/q h eQ h 3 v G H^Q) : divv = q h and ||v|| 1)n < ^\\q h \\o,n- 
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It is well known that for all the families considered here, an interpolant v 1 G Vh of v exists 



that verifies (2.19) (in particular for s = 1), and 

divi/ = divu (= qh). 

By observing that \v\ = on the internal edges as v G H 1 ^), and by using the Agmon 
trace inequality [2] and ( |2.19[ ) (for s — 1), we have 

(3.12) \lv J ]\l := KVtlWl = K'WUv 1 ~ v) t j\\l < C \v\l n . 



Hence, again using (2.19), we deduce that 

||« 7 ||dg < C \v\ i a . 



Thus (3.11) is proved. 



To prove (3.10) we need a discrete Korn inequality, which we prove in the following lemma. 



Also see Appendix [A] for further comments on the validity of the result in three dimensions. 

Lemma 3.2. Let Vh be a piecewise polynomial subspace of Ho^ n (div;Q). Then, 3Ca- > 
independent of h such that 



(3.13) 



\v\l, h <c K IK^fc + ^r 1 ! 



V 



* Jl No,. 



Proof. To show (3.13), a direct application of P, Inequality (1.12)] to v G Vh gives 
(3.14) \v\l h < C 



K 



\e(v) 



\o,T h 



\ 



E /l e 1 H[«t]|lS,e+ SUp (fv-Tjds 

ee£0 rieRM(n) \Jr 

ll T ?llo,r=l, j r r?=0 



/ 



where RM(Q) is the space of rigid motions on Q defined by 

RM(tt) = {a + bx : a G M 2 be so(2) } 
with so(2) denoting the set of skew-symmetric 2x2 matrices. 



We now show that the last term in (3.14) can be bounded by the first two. We claim that 
(3.15) sup (fvrids) <c(\\s(v)\\l Th + Y,K 1 \\lv t \ 



r]eRM(n) 
IMIo,r =1 , frO=0 



IL ; • 



There are surely many ways of checking (3.15). Here, we propose one. For v G Vh and 
r) G RM(£l) with J r r)ds = 0, we set 



X(v, rj) := J v ■ rids, 



DG FOR STOKES 



11 



and we want to prove that 

(3.16) T(v, tj) < C(\\e(v) g i% + £ K 1 ||[^] ||g >e ) 1/2 ||r7||o,r 



will easily give (3.15) taking the supremum with respect to r\ with ||r/||o,r = !• To prove 
(3.16) for every r\ £ RM(Q) with j r rjds = 0, we consider the following auxiliary elasticity 
problem: 





\ T 


= e( X ) 


in Q 


(3.17) i 


divr 


= 


in Q 




r ■ n 


= ~n 


on r. 



The above problem has a unique solution x provided we complement problem (3.17) with 
the side condition 

f x -Cdx = o VCe#M(ft). 

Jn 

Due to well-known results on the regularity of the solutions of PDE systems on polygons, 
the solution r of (3.17) (which, a priori, on a totally general domain would only be in 
{L 2 {Vt)) 2 S y^) satisfies the following a priori estimate: a p > 2 (depending on the geometry of 
Q) and a constant C p exist such that for all rj £ RM(Q) the corresponding r satisfies 



(3.18) 



\ T \\(Lv(Q))m 



< C p ||77||o,r- 



The following technical proposition is the 2-dimensional version of Proposition A.l given 
in Appendix |X} The result stated here is, restricted to the case of interest, d = 2: 

Proposition 3.3. Let T be a triangle with minimum angle 9 > 0, and let e be an edge ofT. 
Then for every p > 2 and for every integer k maX! a constant C Pi g^max exists such that 



(3.19) 



-1/2 — 

v ■ (r • n) ds < C pfijkmax h T \\ v|| , e h T v ||r|| ,p,r 



for every r £ (L p (fi))^^ having zero divergence and for every v £ p kmax (T). 

Integration by parts together with the second equation of (3.17), ( 2.4[ ) and (3.5), gives 



(3.20) 



l(v, r]) = J v ■ rjds = J v ■ (r • n)ds 

= (e(v) ■■ r) Th + (^,divr) Th - {{v t \ : {r}) £ ° 
= {e{v):T) Th -{\v t \:{T}) £ o. 



At this point, we can apply (3.19) to each e of the last term in (3.20) and then use the 
generalized Holder inequality (with q — 1/2 and r = 2p/(p — 2), so that - + - + ^ = 1). 
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Therefore, we obtain 

(3.21) flM--{r}ds< cf^K^v 



i il lo,, 



1/2 



l0,p,T(e) 



l/P 



p— z 

5>- 



l/r 



< c|K]|jT|| 0iPi nMfi) 1/r 



where for each e & with e = <9T + n dT~ , the set T(e) refers to T(e) := T + U T~ . In 
the second line, p(Q) denotes the measure of the domain Q, whereas the constant C still 
depends on p, k max and on the maximum angle in the decomposition Th- 



From (3.20), (3.21), and the bound (3.18) we then obtain 
(3-22) |I(«,»l)l<C(l|e(*)||o,r k + |[«*]U)ll»lllo^ 



which gives (3.16). Thus the proof of the lemma is complete. 



□ 



Remark 3.4. The fact that in inequality (3.13) only the jumps over the interior edges 



e E (but not on the boundary edges) are included, prevents a direct and straightforward 
application of the results from [6] . The proof presented here is surely too elaborate, and we 
believe that a simpler proof is possible. However some of the machinery used here is likely to 
be of use elsewhere. Therefore, we decided that it would be worthwhile to present the proof 
we have obtained to date. 



The stability of dh(-, •) in the 
machinery. We have 



|i3G- n orm can now be easily checked with the usual DG 



{e(v)}:[v t ]ds 



< h l ' 2 \ 



e (v 



0,e 



\h-W[v t 



0,e, 



which when we proceed as in [3] (or as in (3.21 ) with p = 2) yields 



(3.23) 



J2 f{e{v)}:[v t ]da 



<c>UIK]l*- 



Using (3.23) in (3.7), we then have 

a h (v,v) > 2u\\e(v)\\l Th + 2u oc\\[v t ]\l - AuC\v\ llTh \[v t \ 



Now using the Korn inequality ( |3 . 1 3 ) and the usual arithmetic-geometric mean inequality, 
we easily have a big enough a : 

We close this section with the following theorem. 



Theorem 3.5. Let (Vh, Qh) be as in one of our three examples. Then problem (3.6) has a 
unique solution (uh,Ph) € Vh X Qh that verifies 

(3.24) divu h = in il. 
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Moreover, there exists a positive constant C , independent of h, such that for every Vh G Vh 
with diwh = and for every qh G Qh the following estimate holds: 



(3.25) \\u - u h \\ DG < C \\u - v h \\ DG , 



\p - Ph\\o,n <C(\\p- q h \\o,n + \\u - v h \\ DG ), 



with (u,p) solution of (1.5) 



Proof. The existence and uniqueness of the solution of (3.6) follow from (3.10 )-( 3. 1 1 ). The 



divergence-free property (3.24) is implied by (2.14), which holds for all our choices of spaces. 



Let Vh G Vh also be divergence- free; then we obviously have that b{vh — Uh, q) — for 
every q G L 2 (f2)/IR. In particular, b{v h — u h ,p — p h ) = 0. Hence, from the coercivity (3.10), 



consistency (3.8), and continuity of cih(-, •) we deduce immediately 



j\\v h - u h \\ DG < a h (v h - u h , v h - u h ) = a h (v h -u,v h - u h ) < \\v h - u\\ DG \\v h - u h \\ DG . 



On the same basis we deduce that the first estimate in (3.25) follows by triangle inequality. 
For every Wh G Vh, using the consistency and continuity of a^(-, •), we have 

b(w h , q h - p h ) = b(w h , q h -p) + b{w h , p - p h ) = b{w h , q h -p)- a h (u - u h , w h ) 

(3.26) < {\\q h -p|| ,n + \\ u ~ u h \\ DG )\\w h \\ DG . 



By dividing (3.26) by ||iUft||ix? and then using the inf-sup condition (3.11), we immediately 
deduce that 

P\\Qh ~ Ph\\o,a < Uh ~ p\\o,n + \\u - u h \\ DG , 
and that the second estimate in (3.25) follows again by triangle inequality. □ 



Remark 3.6. In the assumptions of Theoren \3.5 , we could obviously consider any trio of 
finite element spaces satisfying HO. However, for choices like RT 0; not considered in our 



three examples, the estimate (3.25) could be meaningless, as the term \\u — Vh\\DG does not, 



in general, go to zero with h. Still, this choice could be profitably used, in some cases, as a 
preconditioner, as it does satisfy HO, HI, and H2. 

4. Discrete Helmholtz decompositions 



So far, we have assumed that the computational domain Q is a polygon. From now on, 
for the sake of simplicity, we are going to work under the stronger assumption that Q is a 
convex polygon. As is well known, this allows the use of better regularity results, and in 
particular the if 2 -regularity for elliptic second-order operators. 

Following [H] we define the discrete gradient operator Q h : Q h — > Vh as 

(4.1) (Ghqh,v h )o,n ■= -(q h ,div v h )o,n Vi^ G Vh- 

Lemma 4.1. Assume that together the three spaces (Vh, Qh,Nh) satisfy assumption HO 



(given in Definition 3.1). Then for any Vh G Vh a unique qh G Qh and a unique (fh G Nh 
exist such that 



(4.2) 



v h = Qhqh + curliph, 
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that is, 

(4.3) V h = g h {Q h ) © curlN h . 

Moreover, a constant C independent of h exists, such that the following estimate holds: 

(4.4) WQhQhWDG < C||div« fc || ,n- 
Proof. For Vh € Vh, consider the auxiliary problem: 

dq 



(4.5) 



Aq = divvh in fi, 



dn 



on <9f2, and / qdx = Q. 



Owing to the boundary conditions in Vh, we have that divvh has zero mean value in Q. 

? solution, that satisfies 

\qh,n < C reg || dwv h || o,n- 



Hence, problem (4.5) has a unique solution, that satisfies 
(4.6) 



We write (4.5) in mixed form: 

(j = —Vq in Q, diver = divt^ in Q, cr • n = on dQ. 

and, we consider directly the approximation of the mixed formulation: Find (<Th, qh) £ 
^ x Q h such that : 



(4.7) 



(fft, T) ,n - (<?h, div r) 0) n = 

(divert, s/Jo.n = (div s ft ) ,n 



(4- 



Problem (4.7) obviously has a unique solution, which moreover satisfies 
||cr - cr h || ,Q < Ch \cr\i,n < CC reg h \\divv h \\ ,n, 



given that (4.6) was used in the last step. As both Vh and cr^ are in Vh (and as (2.14) 
holds,) the second equation in (4.7) directly implies that 

div (a h - v h ) = 0. 



Hence, the exact sequence (|2.8|) implies that 
(4.9) 



a unique (ph £ -A4 exists such that CTh — Vh — cur\(p h - 



Next, by using the first equation in (4.7) and then applying definition (4.1), we deduce that 
((?h,T)o,n = (<?/>, div r) ,n = -{QhQh, f)o,n Vr G V fc , 



which implies <Jh = —Qhqh, that joined to (4.9) gives (4.2). 
In order to prove (4.4), we recall that 

(4-10) \\Ghqh\\DG= \Wh\?DG 



l v<T fcllo,T h 



{<7h)t_ 
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For the first term, by adding and subtracting the interpolant a 1 of cr and then using inverse 
inequality and (2.19), we have: 

||V<r h ||o,r h < l|VK-<7 / )|| 0) r, + IIV^Ho^ 

(4.11) < C inv hr l \\(T h - O-'Ho.Tfc + C||V<T|| ,T h . 



From triangle inequality, (4.8), and standard approximation properties (see (2.19)), we have 
(4.12) \W<T h \\ QJh < C||div«fc|| 0l n. 



The jump term in (4.10) is estimated similarly. First, we remark that cr = — Vg with 
q G H 2 (Q) so that [cr] = on each e G therefore 



Then, using Agmon trace inequalities (|4.8|) and the boundedness of ah and cr, we have 

1 2 



\cr h t - cr, 



{<T h ) t -(T t J||| 



< C t /i- 2 ||cr h - cr || 2 Tft + C t \\V((T h - (T)\\l Th 

< CC reg \\divv h \\l n . 

Thus the proof is complete. □ 

5. Preconditioner: Fictitious Space Lemma and Auxiliary Space Framework 

5.1. Preconditioner for the semi-definite system. Assume V is a Hilbert space equipped 
with the norm || • ||v and that A : V t— > V is a bounded linear operator. We define the bilinear 
form 

(u,v)a = (Au,v). 

We say A is symmetric if the bilinear form (u, v)a is symmetric. We say that A is semi- 
positive definite if 

{v,v) A >0, VveV 

and a > exists such that 

{V, v)a > Ot\\v \\y/jf(A)> 

\/v G V/N{A). 

And, A is SPD if a > exists such that 

Cu>f),4 > oi\\v \\y, Vf G V. 
One useful property of the symmetric a semi-positive definite operator is that 
(5.1) Av = iff (Av, v) = 0. 

A preconditioner for A is another symmetric semi-positive definite operator B : V' i— >■ V. 
Again, we consider the bilinear form 

tf,g) B = (f,Bg). 
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The operator BA : V i— > V satisfies 

(BAu, v) A = (Av,BAu) = (Au,Av) B - 

Lemma 5.1. If A : V H- V and B : V' i— >■ V are both symmetric semi-positive definite such 
that B is positive definite on R(A), then 

(1) B : R(A) i — y R(BA) is an isomorphism (with the inverse satisfying trivially that 
B- 1 (BAv) = Av). 

(2) The bilinear form (-,-)b-i defines an inner product on R(BA). 

(3) The bilinear form {-,-)a defines an inner product on R(BA). 

(4) BA is symmetric positive definite on R(BA) with either of the above two inner prod- 
ucts. 

Proof. All these results are pretty obvious, and their proofs are similar. Let us give the proof 
for 3 as an example. 

We only need to verify that (•, -)a is positive definite on R(BA). If v G R(BA) is such 



that (v , v)a = 0, then, by (5.1 ), we have Av = 0. We write v = BAw for some w E V, then 
ABAw = and hence (Aw, Aw)b = 0. As B is positive definite on R(A), we have Aw = 0. 
Thus, v = ABAw = 0, as desired. □ 

For the system Au = f, we can apply the preconditioner B and the preconditioned con- 
jugate gradient (PCG) method with respect to the inner product (■, with the following 
convergence estimate: 

ii„_ IliA < 2 fv™^yn„_„.| U . 

The condition number can then be estimated by k(BA) < ci/cq, either where 

Coiy,v) B -i < (BAv,v) B -i < Ci(v,v) B -i, \/v e R(BA), 
or equivalently where 

c (w,w)b < (Bw,Bw)a < ci(w,w)b, Vu> G R[A), 

or where 

Ci\v,v)a < {B^v^v) < c^ 1 (v,v) A Vt» G R{BA). 

5.2. Fictitious space lemma and generalizations. Let us present and prove a refined 
version of the Fictitious Space Lemma originally proposed by Nepomnyaschikh [16] (see also 

TO- 



Lemma 5.2. Let V and V be two Hilbert spaces, and let LT : V i— )■ V be a surjective map. 
Let B : V i-> V be a symmetric and positive definite operator. Then B := UBU' is also 
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symmetric and positive definite (here IT : V i— > V"' is swc/i i/m£ (II'p, {;) = (g, IK), g E V 
and v EV). Furthermore, 

(B^v, v) = inf (B~ 1; u, v). 

Tlv=v 

Proof. It is obvious that B is symmetric and positive semi-definite. Note that if v G V" such 
that = 0, then (SII'v, iTu) = = 0. This means that U'v = as B is SPD. 

Hence, u = as IT is injective. This proves that B is positive definite. 

For any v G V, let v — Hv and v* = BTl'B~ l v. As we obviously have Hv* = v, we can 
write v = v* + w with Hw = 0. Thus, 

inf (B~ x v,v) = inf (JET 1 ^* + w), v* + w) 
nii=u riui=o 

= (fi- 1 i)*,i5*) + inf ((B^w^w) +2{B- 1 v*,w)) 

From the definition of v* we have 

(B- 1 ^^*) = (B-\Uv*) = (S _1 t;,t;), 

and also 

(B^v^w) = (B^BU'B^v^w) = (WB- X v,w) = {B~\mb) = 0. 
The last two identities lead to the desired result. □ 

Theorem 5.3. Assume that A : V h- > V' and A : V i— >■ V' are symmetric semi-definite 
operators. We assume that U : V ^ V is surjective and that U(N(A)) = N(A). Then for 
any SPD operator B : V' i— > V, we have, for B = UBU' , 

k{BA) < k(U)k{BA). 
Here k(IT) is the smallest ratio c\/cq that satisfies 

(5.2) c^{Av,v) < inf (Av,v) < c^ l {Av,v), \fv G R{BA). 

Hv=v 

Proof. Denote k(BA) = 61/60 with 61 and 60 satisfying 

b^(v,v) A < (B- l v,v) < b \d,v) A , W G R{BA). 



By (5.2), we obtain 



6r 1 cr 1 |MlA < inf (B^v,v) < 6 1 c 1 |H|^, Vt> G R(BA). 

Uv=v,v€R(BA) 

By the assumption that U(N(A)) = N(A), we can prove that W(R(A)) C R(A) and 
{v\Uv = ve R(BA)} = {v\Uv = ve R(BA),v G #(£?!)}. 



By Lemma 5.2 



inf (B^v, v) = inf (B'^v) = (B~\v), Vt> G R(BA). 

Tlv=v,v£R(BA) llv=v 



18 

Therefore, 



B. AYUSO DE DIOS, F. BREZZI, L. D. MARINI, J. XU, AND L. ZIKATANOV 



&rVNft < (B-\v) < &bVNft Vt; G fl(SA). 



□ 



Theorem 5.4. ^sst/me t/iat the following two conditions are satisfied for IT. First, 

\\Rv\\a < c i||^||a) ^ e ^- 
Second, for any v £ V , v £ V exists such that Hv = v and 

\\v\\a < c olMU- 



Then k(U) < Ci/cq and, under the assumptions of Theorem 5.3, 



k{BA) < 



k{BA). 



Remark 5.5. In view of the application of the above results to our case (as we shall see 
in the next subsection), it would have been enough to restrict ourselves to the symmetric 
positive definite case (instead of the semi-definite case treated in the last two subsections). 
However we preferred to have them in the present more general setting, as in this form they 
are likely to be useful in many other circumstances (starting, as natural, from the extension 
of the present theory to the three-dimensional case). 

5.3. Application to our problem. In this section we design a simple preconditioner for 
the linear system resulting from the approximation of the Stokes problem (1.5) defined in 
(3.6)-(3.7). Note that the bilinear form ah(-, -) defined in (3.7) provides a discretization of 
the vector Laplacian problem 

— d\v{2u£{u)) = f in f2, u ■ n = 0, (s(u) ■ n) • t — on T. 



We denote by A the operator associated with a^-, •). As the solution u h e V \ of (3.6) is 
divergence-free, the discrete Helmholtz decomposition (4.2) implies that 

a unique iph G Afh exists such that u h = curl-?/v 

o 

At this point, it is convenient to introduce the space Vh as 
(5.3) V h :=V h nH (dw°;Q). 



And, we note that as the sequence (2.8) is exact, we have 
(5.4) V h = curl A/;, 

and that the mapping is one-to-one. Therefore, restricting the bilinear form a^( 
in the spirit of Remark 2.1 



to V h 



corresponds here to restricting the trial and test space to V \= 
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curl(A//i). The discrete problem (3.6) then reduces to the following problem: Find iph GV/j 
such that 

o 

(5.5) a h (tp h ,cp h ) = {f,(fh) V<p h eV h 

o o 1 o 

Defining the operator A h :V h ^V h by (A h ip h , (p h ) = a h (ip h , (p h ), ip h , ip h eV h , we can write 
(5.5) as 

A h ip h = fh- 

o _ 

We now use the original space V \ as the auxiliary space for Vh- Define A h : V \ i— > V' h by 

(A h Uh,Vh) = a h (uh,Vh),Uh,Vh G Vh- We note that A h is a discrete Laplacian. We assume 

that Bh is an optimal preconditioner for A^. 
We now define the operator 



(5.6) 



cur\{N h ) 



which maps each v\ t G Vh to its divergence free part in Vh, according to (4.2). Note that 11^ 

o 

is a surjective operator and that n fe acts as the identity on the subspace Vh- The auxiliary 
space preconditioner for Ah is then defined by 

(5.7) B h = U h B h n* h . 

Lemma 5.6. Assume that the spaces(Vh, Qh,Nh) satisfy assumption HO. Then Bh given 



by (5.7) is an optimal preconditioner for Ah as long as Bh is an optimal preconditioner for 
A h - 



Proof. Following the auxiliary space techniques (Theorem 5.4), we need to check that the 
following two properties are satisfied: 

(Al): Local Stability: there exists a positive constant C\ independent of h such that 

(5.8) Hn^/Jix; 

< Ci\\v h \\ dg e Vh 

(A2): Stable decomposition: there exists a positive constant C2 independent of h such 

o 

that for any Wh GV/j there exists v h G Vh such that Iih v h = w h and 

(5.9) < ^II^Udg- 



To prove (5.8) from the Helmholtz decomposition (4.2) and the definition (5.6) of n^, we 
have 



(5.10) 



v h = QhQh + cuvltph = QhQh + HhVh- 



Using estimate (4.4) from Lemma 4.1 and the clear fact that divu^ is the trace of e(vh), we 
have 



(5.11) 



\Gh.qh\\DG < C||divu h || 0l n < C||e(t;)||o,r h < C\\v h \\DG- 
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Hence, (5.8) follows from (5.10) and (5.11): 

lln^Hzx? = \\Vh - OhQuWdG < \\Vh\\DG + WQhQhWDG < C\\v h \\ DG . 



Finally, the inequality (5.9) holds with = 1 by taking = Wf 



□ 



6. Numerical experiments 



6.1. Setup. The tests presented in this section use discretization by the lowest order, 
namely, BDM4 elements paired with piece-wise constant space for the pressure. They verify 



the a priori estimates given in Theorem 3.5 and confirm the uniform bound on the condition 



number of the preconditioned system for the velocity. 

As previously set up, the discrete problem under consideration is given by equation (3.6) 
with bilinear forms a/ l (-, ■) and &(•, ■) defined in (3.7). In the numerical tests presented here, 
we take v — 1/2 and the penalty parameter a = 4 in (3.7). We present two sets of tests 





(a) Coarsest mesh 



(b) Mesh for level of 
refinement J = 3 



Figure 6.1. Meshes used in the tests for the unit square domain 

n = (0,1) x (0,1) 

with A corresponding to the Stokes equation discretized on a sequence of successively refined 



unstructured meshes as shown in Figures 6.1-6.2 On the square the coarsest mesh (level 
of refinement J = 0) has 160 elements and 97 vertices with 448 BDM degrees of freedom. 
The finer triangulations of the square domain are obtained via 1, . . . , 5 regular refinements 
(every element divided in 4) and the finest one is with 163, 840 elements, 82, 433 vertices 
and 490, 496 BDM degrees of freedom. Similarly for the L-shaped domain we start with 
a coarsest grid (J = 0) with 64 vertices and 97 elements. For the L-shaped domain the 
finest grid (for J = 5) has 99, 328 elements, 50, 129 vertices and 297, 056 BDMx degrees of 
freedom. In the computations, we approximate the velocity component Uh of the solution of 
the Stokes equation by solving several simpler equations (such as scalar Laplace equations). 
After we obtain the velocity, the pressure then is found via a postprocessing step at low 
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(a) Coarsest mesh (b) Mesh for level of 

refinement J = 3 



Figure 6.2. Meshes used in the tests for the L-shaped domain 

n = ((o J i)x(o,i))\([| ) i)x[i i)) 



computational cost. Further, for this sequence of grids the BDMi interpolant of a function 
v on the k — th grid is denoted by v Ik . Accordingly the piece-wise constant, /^-orthogonal 



projection of p is denoted by p Ik . We also use the notation (uk,Pk) for the solution of (3.6) 
on the k — th grid, k = 0, . . . , 5. 

6.2. Discretization error. We now present several tests related to the error estimates 
given in the previous sections. We computed and tabulated approximations of the order of 
convergence of the discrete solution in different norms. These approximations are denoted 
by 7o ~ Po, Jdg ~ Pdg, lv ~ Pv an d 7* ~ P*- The actual orders of convergence /3 , Pdg, 
P p , and p* are 

\\u - « h || 0j n ~ C(u)h^\ \\u - u h \\ DG w C(u)h^, 
||p-Ph||o,n « C(u,p)h^, \[u h \U « tf(u)/^. 



Here, as in (3.12), we denote 



Note that /3* is the order with which the jumps in the approximate solution (not in the error) 
go to zero. 



We present two sets of experiments to illustrate the results given in Theorem 3.5 First, 
we consider the exact given solution and calculate the right-hand side and the boundary 
conditions from this solution. We set 



(6.1) </> = xy(l-x)(2x-l)(y-l)(2y-l), u = curk/>. 
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Square domain 


k 


1 


2 


3 


4 


5 


7o 


1.56 


1.79 


1.91 


1.96 


1.98 


IDG 


1.32 


1.58 


1.66 


1.66 


1.62 


% 


0.43 


0.70 


0.86 


0.94 


0.97 


7* 


0.70 


0.85 


0.94 


0.97 


0.99 



L-shaped domain 


k 


1 


2 


3 


4 


5 


7o 


1.43 


1.64 


1.84 


1.93 


1.97 


IDG 


1.17 


1.47 


1.60 


1.62 


1.60 




0.24 


0.52 


0.79 


0.91 


0.96 


7* 


0.64 


0.76 


0.89 


0.95 


0.98 



Table 6.1. Approximate order of convergence for the difference (u 1 — Uh) and 



(p 1 — p h ) and the jumps 



u and p are given in (6.1) and (6.2) 



UhjL for the square and L-shaped domains. Here, 



Clearly, the function <p vanishes on the boundary of both the domains under consideration 



and we take u defined in (6.1) as exact solution for the velocity for both the square and the 
L-shaped domains. For the pressure we choose as exact solutions functions with zero mean 
value and select p different for the square and the L-shaped domain, namely 



(6.2) 1 



p = x 2 — 3y 2 + (square domain), 

p = x 2 — 3y 2 + yxy, (L-shaped domain). 



The right hand side / is calculated by plugging (u,p) defined in ( |6. 1 )— ( 6.2 ) in (1.1 ). Table 6.1 



shows tabulation of the order of convergence of (uh,Ph) to (u 1 ,p ! ) for both the square domain 
and the L-shaped domain. The values approximating the order of convergence displayed in 
Table 16.11 are 

7 = log 2 — n— r ii— > 7* = log 2 



U k - U k\\ \l u k-l\\* 

\\Pk-l — Pft-l||o,fi , 1 r 

l|Pfc-P*llo,n 

Here || • || stands for any of the DG or L 2 norms. The quantity 7 is the corresponding 70 or 
7dg- From the results in this table, we can conclude that in the || • \\dq norm the dominating 
error is the interpolation error, and as the next example shows, in general, the order of 
convergence in || • \\ DG is 1. 

The second test is for a fixed right hand side / = 2(1, x). We calculate approximations to 
the order of convergence of the numerical solutions on successively refined grids as follows: 

\\u k -Ufc_i|| |[«fcll* - |[«fc-l]|* 
7 = log 2 ^ rr, 7*=log 2n in nf IT' 



7 P = log. 



|«jfc+l -ttfell IL«fc+lJII* - \l u k. 

\\Pk -Pfc-i||o,n 
Ibfc+i -Pfc||o,n : 
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Square domain 


k 


1 


2 


3 


4 


7o 


1.77 


1.86 


1.93 


1.97 


IDG 


0.87 


0.92 


0.97 


0.98 


% 


0.98 


0.99 


0.99 


1.0 


7* 


0.90 


0.95 


0.98 


0.99 



L-shaped domain 


k 


1 


2 


3 


4 


7o 


1.28 


1.44 


1.56 


1.50 


IDG 


0.53 


0.46 


0.44 


0.40 




0.97 


0.93 


0.87 


0.75 


7* 


0.40 


0.40 


0.37 


0.35 



Table 6.2. Approximate order of convergence of the error for square and L- 
shaped domains and right-hand side / = 2(1, x). 



Again, || ■ || denotes any of the (semi)-norms of interest and 7 approximates the corresponding 



order of convergence. Table 6.2 shows the tabulated values of 70, 7dg ; 7p> an d 7*. It is clear 
from these values that the order of approximation for the velocity and the pressure is optimal 
for the square domain, whereas for the L-shaped domain the convergence is not of optimal 
order, due to the singularity of the solution near the reentrant corner. 

6.3. Uniform preconditioning. The tests presented in this subsection illustrate the effi- 



cient solution of the system (6.3) below by Preconditioned Conjugate Gradient (PCG) with 



the preconditioner given in (6.4). We introduce the matrices representing the bilinear forms 



defined in (3.6)-(3.7), and also the mass matrix for the BDMx space. We denote by M the 



mass matrix on V h and by A the stiffness matrix associated with %(•, •) on Vh in (3.6) 



(3.7). We note that A, without the divergence-free constraint, is spectrally equivalent to 



in (3.6) is made of vector fields that are curls of 



two scalar Laplacians. 

It is known that the null space of b(- 
continuous, piecewise quadratic functions vanishing on the boundary. We denote by P cur i 
the matrix representation of the natural inclusion of these curls in the BDM space. Namely, 

curl(basis functions in Mh) = (basis functions in Vft)P cur i. 

It is easy to see that 

Ag = P curl MP cur l. 

where A q is the discretization of the Laplacian on with homogeneous Dirichlet boundary 
conditions. 



The problem of finding the solution of (5.5 ) then amounts to solving the following algebraic 
system of equations 



(6.3) 



■PcurlAPcurlTJ 



curl ■ 



Here the superscript T means that the adjoint is taken with respect to the £ 2 -inner product, 
U is the vector containing the velocity degrees of freedom, and F is the vector representing 



the right-hand side (f,v) of the problem (3.6). 
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The matrix representation B of the preconditioner B described in the previous section has 
the following form: 



(6.4) 



B 



A-^MBMPchA- 



where B is a preconditioner of A. In summary, the action of the preconditioner requires 
the solution of systems corresponding to 4 scalar Laplacians. It is also worth noting that 
suitable multigrid packages for performing these tasks are available today. 



The results are summarized in Table 6.3 The legend for the symbols used in the table 
is as follows: n it is the number of PCG iterations; p is the average reduction per one such 



iteration defined as p 



\\ r n u \\l2 



; J is the refinement level, for which h « 2 J ho, where 



ho is the characteristic mesh size on the coarsest grid. From the results in Table 6.3 



we 



Square domain 


J 





1 


2 


3 


4 


5 


n it 


5 


6 


6 


7 


7 


7 


P 


0.06 


0.07 


0.09 


0.10 


0.11 


0.13 



L- shaped domain 


J 





1 


2 


3 


4 


5 


n it 


6 


7 


8 


9 


9 


10 


P 


0.08 


0.10 


0.13 


0.18 


0.21 


0.21 



Table 6.3. Preconditioning results for square domain (left) and L-shaped do- 
main (right). The PCG iterations are terminated when the relative residual is 
smaller than 10 -6 . 



can conclude that the preconditioner is uniform with respect to the mesh size. It is also 
evident that this method is in fact quite efficient in terms of the number of iterations and 
the reduction factor. 

To conclude this section, let us point out that the generalizations to the 3D version of the 



preconditioner, and the introduction of the approximate solvers in (6.4) can be achieved in 
a similar fashion and are the subject of ongoing research. 



Appendix A. Proof of Proposition 13.31 
We now state and prove a result, Proposition [AT given below, which contains as a partic- 



ular case the Proposition 3.3, used in Section [3] to show Korn inequality (cf. Lemma 3.2) for 
the case d = 2. After giving its proof, we comment briefly on how the result can be applied 
to show the corresponding Korn inequality (3.13) (cf. Lemma 3.2) for d = 3. 



Proposition A.l. Let T be a triangle (or a tetrahedron for d = 3) with minimum angle 
9 > 0, and let e be an edge (resp. face) of T. Then for every p > 2 and for every integer 

fcmax Q> Constant Cpft^kmax 

exists such that 
(A.l) / v ■ (r • n) ds < C pAK 



h T IM|o,e h T 2P 



\ T \\0,p,T 
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for every r G {L p {VL)) d s ^ having zero divergence and for every v G p kmax (T) 



I sym 



Proof. First we go to the reference element t: 



(A.2) 



v 



(r • n) ds < Cg\e\ 



v ■ (t - n) ds 



<Cah 



d-X 



V 



if ■ n) ds 



where v and r are the usual covariant and contravariant images of v and r, respectively. 
And, here and throughout his proof, the constants Cq and Cg t k max may assume different 
values at different occurrences. Note that the divergence-free property of r will be preserved 
by the contravariant mapping, and that v will still be a vector-valued polynomial of degree 
< k max . Then for every component v of v, we construct the auxiliary function ip v as follows. 

First we define tp v on dT by setting it as equal to v on e and zero on the rest of ST. Then 
we define p v in the interior using the harmonic extension. It is clear that (p v will belong to 

W l,p ' (T) (remember that p > 2 so that its conjugate index p' will be smaller than 2). Using 
the fact that v is a polynomial of degree < k max , it is not difficult to see that 

( A - 3 ) IIVvllwLf'CT) - Co,k m a X \\v\\o,e- 

A simple integration by parts (remembering that div r = 0) now gives 



(A.4) 



v ■ (t ■ h) ds 



/ ip v ■ (t • n) ds = / Vip v : f dx 

J&t Jf 



- l^lw 1 'P'(f)ll''"ll(LP(f))^ - Ce,k max \\v || 0,e 



l(LP(T))! 



The proof then follows immediately by inserting (A.4) into (A.2) and then coming back to 
the current element T by means of: 



(A.5) 



Hide < C e h t 



d-l 
2 



l T ll(Lf(f))f^ - C '^T P |l r ll(LP(T))^ 



□ 



With this result in hand, we can show the Korn inequality (3.13) given in Lemma 3.2 for 
d = 3. It is necessary to modify the proof in only two places: the definition of the space of 



rigid motions on f2, J2Af(0), and the application of Proposition 3.19 The space RM(tt) 
is now defined by: 

RM(tt) = {a + bx : a el rf be so(d) } 
with so(d) denoting the Lie Algebra of the skew-symmetric d x d matrices. 



To prove (3.16) (and so conclude the proof of (3.13)), estimate (3.21) is replaced by 



estimate (A. 6) below, which is obtained as follows: first, by applying (A.l) (instead of 



(3.19)) from Proposition A.l to each e in the last term in (3.20) and then by using the 
generalized Holder inequality with the same exponents as for d — 2 (with q = 1/2 and 
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r = 2p/(p-2),so thati + i + i = l) 

1/2 



(A.6) J2 M: {r}ds< C^^lMl) (E 



l(),p,T(e) 



<p-2) r N 1/r 



< C f |[«t]|2l|r||o^A*(n) 1/r 



Here, as in estimate (3.21 ), T(e) in the first line denotes the set of elements T whose boundary 
contains e, /x(fi) in the second line denotes the measure of the domain Q, and the constant 
C still depends on p, k max , and on the maximum angle in the decomposition T h . The rest of 
the proof of Lemma [372] proceeds as for d = 2. 
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